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Abstract 



■ We present elliptic solutions to the background equations describing the Lemaitre-Tolman- 

W)' Bondi (LTB) metric as well as the homogeneous Friedmann equation, in the presence of dust, 

curvature and a cosmological constant A. For none of the presented solutions any numerical 
• integration has to be performed. All presented solutions are given for expanding and collapsing 

^ I phases, preserving continuity in time and radius; both radial and angular scale functions are 

, given. Hence, these solutions describe the complete spacetime of a collapsing spherical object 

00 ' in an expanding universe, as well as those of ever expanding objects. In the appendix we 

, present for completeness a solution of the Friedmann equation in the additional presence of 

radiation, only valid for the Robertson- Walker metric. 

1 Introduction 

Cosmic structures today have entered the non-linear regime. They can not on all scales be 
described by a linear perturbation theory on top of the Friedmann-Lemaitre-Robertson- Walker 
X ■ (FLRW) metric The simplest step beyond linear perturbation theory is to look at separate 

patches, describing the evolution in each patch as a closed system, insensitive to other perturba- 
tions outside the patch. Over-densities can be studied in the approximation of spherical collapse, 
where underdensities expand and potentially become spherical voids. Either way, these spherically 
symmetric configurations, whether matching to a surrounding FLRW metric or not, arc described 
by the Lemaitrc-Tolman-Bondi (LTB) metric [l-Q' 

ds^ = -dt^ + S^{r,t)dr^ + R^ir,t){dd^ + si-a^ 6 d(j>^), (1) 

with 

Sir,t)^^Eh^. (2) 

The LTB metric reduces to the FLRW metric if one sets R{r,t) = r a{t) and E{r) = — ^kr^. 

In Ref. it was shown that when peculiar velocities are small, a seemingly non-linear solution 
to the metric becomes a linear perturbation on the FLRW metric in the Newtonian gauge. Here 
we focus on solutions to the metric in full generality. 

Spherical collapse is studied for example for the formation of black holes as well as for determin- 
ing, in cosmology, if and when an initially linear over-density produced during inflation collapses 
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and decouples from the background expansion. In simplest approximation, one considers a homo- 
geneous overdosed patch that expands and collapses, matched to an expanding background by a 
singular shell [loj . Choosing a continuous curvature profile in Eq. ^ , allows for an exact solution 
without singular shells. Spherical collapse in either the approximation or the exact approach, 
gives insight in clustering of matter, and thereby has been related to presence of for example Dark 



Energy, amongst other possibilities |lll - l24j . One should note that initial velocities could be such 
that an overdensity evolves to become under-dense, but such a decaying mode corresponds to an 
inhomogeneous Big Bang, which is at tension with the inflationary paradigm in today's favoured 
model of cosmology. 

The formation of voids is studied for two reasons. One reason is the role of voids in the process 



of structure formation j25l - l28{ , the other is the effect that unusually deep under-densities can have 
on our perception of Dark Energy. Some studies consider one large local void (size varying from 
tens to thousands of megaparsecs) 2S-54|, where others consider a distribution of many voids, 
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the so called Swiss- Cheese universe 

In presence of only dust and curvature, analytical solutions to the LTB-equation (which we 
write down later) for t(R, r), i.e. time as a function of local angular scale factor R and coordinate 
radius r. are known in terms of hyperbolic functions (which become trigonometric functions in case 
of complex arguments). However, observations of distant Supernovae, of the Baryonic Acoustic 
Oscillations and of the distance to the surface of last scattering of Cosmic Microwave Background 
photons combined with locally observed expansion rate, demand the presence of a cosmological 
constant [58-61|. It is crucial to realize, although it is not the topic of this paper, that the presence 



of the cosmological constant is only then necessary to explain the observed geometrical distances, 
if one assumes that on large enough scales the universe is still properly described by the FLRW 
metric and one assumes that the angular diameter distance-redshift relation is correctly described 
by the background dynamics only. 

In the presence of a cosmological constant, dust and curavture, the solution for t{R, r) is 
an elliptic integral. Therefore, in most works, authors elude numerical integration of the Ein- 



stein equations, although see Refs. [62,l63|. However, if one for example wants to solve geodesic 
equations, one typically would perform a numerical integration of the geodesic equations over a 
numerical solution of the background. The unknown error in the numerical background solution 
can propagate into the solution to the geodesic equations, possibly leading to unreliable answers, 
or very slow and sometimes unstable codes. 

In absence of a cosmological constant, but in presence of dust and curvature, the solution for 
t{R, r) with positive K{r) is. 
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while for K(r) < 0, the solution becomes a trigonometric in stead of hyperbolic function if one 
propagates the sign of nir) correctly. In this case, one obtains the inverse solution R{r, t) by 
numerically inverting t{R,r), which can be done quickly and at ultimate accuracy, since is 
known by the definition of the Einstein equations. The thus obtained solution is accurate, fast, 
and allows for reliable and fast integration of geodesic equations [H^ . 

The purpose of this paper is to provide all solutions to the background equations t{R, r) that 
have an initial singularity (Big Bang), in presence of dust, curvature and a cosmological constant, 
in terms of Elliptic Integrals in Carlson's symmetric form, which can be numerically evaluated as 
accurate and fast as any elementary function. Compared to the known exact solution in the case 
of no cosmological constant, the solutions presented in this paper are exact, since one eventually 
obtains R{r,t) by quick and reliable numerical inversion of t{R,r). Hence, throughout this paper 
no numerical integration is performed, and solutions are exact but only semi-analytical. See 
for comparison elliptic solutions involving the Weierstrass elliptic function in Refs. and 
references therein. Note that these references only give t{R, r), which in the FLRW case is enough 
to solve for a{t), but which in the LTB case does not suffice: t{R,t) only straightforwardly leads 
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to the angular scale factor R{r,t), but we also present for the first time the radial scale factor 
5(r, <), which is more involved as we shall see later. 

The main improvement in this work in comparison with the existing literature, is the fact that 
we provide solutions for all functions appearing in the LTB metric, involving spatial derivatives 
drR{r^t), necessary for solving for example geodesic equations. We list the solutions for the limits 
where the local expansion transits to collapse, while a neighbouring shell continues to experience 
expansion, at the same time preserving continuity of all functions. Moreover, in the appendix we 
present linear expansions when either the cosmological constant or curvature is small, or both are 
small. 

The solutions presented in this work allow for a plethora of applications. Let us list but a few. 
For example, the solutions can be used for: 

• any numerical work involving a general LCDM background expansion, 

• obtaining the exact metric around a collapsing structure 

• simulating the universe as observed by an observer in a large and deep void, in presence of a 
cosmological constant, allowing for a direct face off between the cosmological constant and 
the void, 

• studying the evolution of voids in structure formation, in a ACDM universe. 

We release a numerical module, written in Fortran, that computes exact metric functions 
and derivatives for a given curvature profile, that can be easily implemented in any code, at 



http://web.physik.rwth-aachen.de/download/valkenburg/ColLambda/ A brief example of 
how to invoke this module is given in Appendix IdI 

This work is organized as follows. In Section [2] we first list for reference the used Elliptic 
Integrals in Carlson's symmetric form. In Section [3] the LTB metric and its Einstein equations are 
discussed. Then in Section|4]we present one of the main results of this paper, being t{a) in terms of 
Carlson's elliptic integrals. Next, in Section[5]we solve analytically for the functions in the metric 
as a function of time t and scale factor a. In Section [6] we provide an example of an application 
of the solutions presented in this work. Finally we conclude in Section [T] In Appendix |X] we 
provide the asymptotic expansions of the solutions and in Appendix |B] we show the solution of 
the Friedmann equation in presence of radiation, matter, curvature and a cosmological constant 
in terms of Carlson's elliptic integrals. 

Throughout this work, square roots of real quantities are taken to be positive, and for all 
fractional powers of complex numbers x we take the principle value of exp(lna;). Extra minus 
signs due to the possible crossing of branch cuts are written explicitly. We use units in which 
Gat = c = 1. Overdots denote time derivatives, primes denote radial derivatives. Our notation 
follows mostly the notation used in for example Refs. 
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2 Carlson's symmetric form of Elliptic Integrals 



Before discussing solutions to the LTB equation, let us list some definitions for completeness. We 
take the following definitions of Carlson's symmetric form of elliptic integrals from Ref. [t^ , 

s{t) ^VtTxy/tTyVtTz. (3) 

RF{x,y,z) =- / — , (4) 



2 7o <ty 

Rc{x,y) =RF{x,y,y) , (5) 
X 3 Z""" dt 

Rj{x,y,z,p) ^- — -, (6) 



2 7o s(t)(t + p)' 

3 f°° dt 
Rn{x,y,z)^Rj{x,y,z,z)^-j^ 

3{xyzy^ =RD{x,y,z) + Roiz.x.y) + (8) 

which are defined for {x, y, z} G C | {x, y, z} ^ (— oo, 0]. These can be evaluated using an iterative 
procedure, up to unlimited accuracy in very few steps, as explained in Ref. [7l| . The definitions 
are valid for complex arguments, and in all these cases at most one argument is allowed to be zero. 



3 Robertson- Walker and Lemaitre-Tolman-Bondi metrics 



The Einstein equations for the FLRW metric and the LTB metric can be written in the same form, 
with the difference that in the former case no other coordinate dependence than time dependence 
is present, and in the latter case the curvature and scale factor are both radius and time dependent 
and radiation is absent. The FLRW metric is, 



ds^ = -dt'^ + a{t)da 



and the Friedmann equation for the FLRW metric is 



nr. 



aHt) a^it) a^{t) 



(9) 



(10) 



where as usual we define today hy t = to, ag = a (to) = 1, Hq = H{tQ). The different components 
and their relative abundances are radiation fir, dust il„j, curvature ilk and the cosmological 
constant ^a. 

The LTB metric is given by 



where 



-dt^ + (r, t)dr'^ + i?2 (r, t) {dO^ + sin^ ( 

^l + 2r2K(r)Af2 
R{r,t) =ra{r,t). 



(11) 



(12) 



(13) 



where is an arbitrary parameter defining the length and mass scales, combined with the choice 
of units Gn = c ~ 1. The Einstein equation leads to the LTB equation 



5 H^l -/f't".')- 



2L{r) 2r^K{r)M 



/r2 



i?3(r,t) i?2(r,t) 



A 

+ 3' 
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where L{r) = dr M' [r)^ \ + 2r'^K{r)M'^ with A/(r) = Att dr S{r,t)R'^{r,t)p{r,t) being the 
total mas|3 inside radius r and p(r, t) is the local matter density. Now there are three functions 
of r that specify the problem, M(r), K{r) and tBB{r), where the latter is the radially dependent 
Big Bang time, tBB{f) = t{a = 0,r). One of these three can be fixed to an arbitrary function by 
redefining the coordinate r — > f = /(r) for some monotonic function /(r), without changing the 
physical description. As discussed in Ref. [lH , none of the three possible coordinate gauges where 
one of the functions is fixed to an arbitrary monotonic function captures all possible configurations. 
In this work we choose the gauge as follows. Demanding a strictly positive matter density, we 
have L'{r) > 0. We choose L(r) = AnM'^r^/S, such that L'(r) > but L'(r) ^ 0. This implies 
that are are no vacuum regions. Then we have 



M'{r) 



1 + 2r'^K{r)AP 



R'{r,t)R^{r,t) 
In this coordinate gauge the LTB equation becomes. 



1 



3K(r) 



A 



(14) 
(15) 

(16) 



and the configuration is completely specified by the two functions nir) and tBsiT). The short- 
coming of this gauge is, as mentioned above, that it does not allow for solutions with true vacuum 
over a non-zero range in r, for which M'(r) = such that K,{r) — > oo. 

3.1 Normalization 

Normalizing a(r*,to) = 1 at a chosen {r*,to}, we have 

47r r!,(n) 



Kb 



3 1 - ^K.{r*) - 17a (n) 



^^2^3iJ^('^*,^o)-A 1 



to EEi(a(r*,to)), 



3Kb 

4-n- 



(17) 

(18) 
(19) 



which is regular for A — > and fJfc — > 0. Also, we write H^, = H^{r^,,to). One can choose an 
arbitrary at which to normalize, but one has to fix it once and for all. Clearly, in the FLRW 
case, the choice of is irrelevant since the r-depence of H(r,t) vanishes, and we find Hq = H^. 



3.2 Towards solving for t{a, r) 
One can define 

fiA(r) 



3i/2(r,io) 
2AP 



K{r), 



1 A 



(20) 
(21) 
(22) 



^ L{r) is the active gravitational mass, while M{r) is the total rest mass. However, in cosmologically relevant 
scenarios 2r'^K{r)hP < 1, such that L{r) ~ M{r). 
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in which case one retrieves the Friedmann equation in absence of radiation when one drops the 
r-dcpendcnccs in equation (|16|) and when one normalizes a(r, to) = a(^o) = 1- For a- generic matter 
distribution in the LTB metric one however has a{r, to) ^ 1. Then, in the LTB metric, these three 
quantities are the relative content at a{t,r) = 1 in a shell at a given radius: dust, curvature and 
cosmological constant respectively. At to the relative matter quantities are then depending on the 
value of a(r, to), e.g. rtmO-~^{r,to) for matter. The reader should note that in this sign convention 

• rjfc > corresponds to an open universe 

• and vice versa < corresponds to a closed universe. 

From this point on we will neglect radiation, although the reader is referred to Appendix |B] 
for a discussion of the solution in presence of radiation, only valid in the FLRW metric. 

Writing A = a{r, t), the general solution to t for the metric (both LTB and FLRW) is given by 
the integral, 

where we take the positive square root and the integral is performed at constant r. Note that the 
Big Bang time tBsif) acts as an integration constant for the left hand side of this equation. 

In the case of existence of at least one positive real root of the polynomial in the denominator 
{i.e. an over-closed universe (FLRW) or over-closed shell at radius r (LTB)), there are two solutions 
for t(A), one for the expanding and one for the collapsing phase. For example labeling the smallest 
positive real root U{r), such that the turning point in the expansion history lies at a{r,t) = U{r), 
the second (collapsing) solution is 

H4t{A)-tBB{r)] = 



where the sign in front of the second integral is determined by whether a changes sign (collapse) 
or not (continued expansion) at a = U . Throughout the rest of this work, we discard cases where 
a does not change sign, i.e., we only consider cases where a is non-zero at a = U{r). The case 
with more than one positive real root then becomes irrelevant: when a changes sign, by symmetry 
the contraction is identical to the expansion, and we only consider the branch of solutions that 
experiences an initial singularity (Big Bang). As a shrinks again towards 0, for each value of a 
the time derivative is exactly —1 times the time derivative for the same value of a during the 
expansion. The higher roots are never encountered. 

4 Solving for the expansion rate 

We rewrite Eq. ^ to 

[t{A) - tBB{r)] / (25) 

for n = 3, where we rewrote the polynomial 

^m{r) + ^K{'r)a + VLf^{r)a^ = ^(^{a - yi){a - y2){a - j/s), (26) 
with yi the solutions to 

S2A(r) S2A(r) 
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Wc now turn Eq. (j25|) into Carlson's symmetric form by subsequently making the change of 
variables a — > c = ^ and c — > & = c — 4- , such that we get 



H,[t{A)-tBB{r)] 



1 a— da 



1 r c-'-^-^dc 



1 (-1)-^ f°° dc 



1 i-l)-^ f°° db 



(28) 



These transformations are valid, since no rotations in the complex plane are involved, and none of 
the roots is transformed onto the path of integration; no branch points come to lie on the positive 
real axis. At most one root will be at the zero of the integration domain in the variable 6, when 
we integrate /j'™ da. This exception is allowed in the definition of the symmetric elliptic integrals 
in Section O 

In terms of physics this is straight forward to see: for any real negative y,„, we have ^ ~ ^ > 0, 
i.e. the branch point lies on the negative real axis of b; for any real positive j/m, we also have 
— ^ > 0, since we integrate at most up to A = The scale factor never grows beyond its 
smallest maximum for t G M. 

We now obtained the solution for the time as a function of scale factor expressed as a symmetric 
elliptic integral in Eq. (pS)) . that is, 

r , ,x , M 2 (-1)"^ f I 11 11 1 1 \ 



3\/ni ffTs ~ ■ \A yi' A y2 A ' ^ 

Y llr)i=l tfm 

One of the first three arguments of Rj{x^y^ z^p) is allowed to be zero, such that the limit of 
A ^ yi for positive real yi is trivial. 

In order to keep a connection to the well-known Friedmann equation, we so far used a notation 
in terms of iJ* and il,i(r). However, as iJ* = H{r,tQ) is a function of r, it is more convenient to 
go back to the original form of Eq. and recast Eq. (^5)) as 

1 /"^ \/a,dd 



t{A)-tBB{r)=^ (30) 



^-^ 'T-— 'T ' (^1) 



VSA /ffi ~ zi A Z2 A z^' A 

V llm=l 



where the roots yi and Zi are related by a simple rescaling. 

4.1 The roots Zi 

The roots Zi are trivial. 



^1 =rT^ - ^iTT-, (32) 

2333Z (f>3 

_^1W3^ _ (1-^V3)$^ 

2i33$l 2l3iz ' ^ ^ 

.3^^i^^-(ii^, (34) 
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where $ = ^/3^/27X2Z4 + 4^3 ^3 _ gxZ^, and X ^, F = 2K(r) and Z = For the FLRW 

metric the roots yi obey the same expressions, replacing X, Y and Z by the corresponding fi^. 



5 The metric functions and their time derivatives 



In this section we present the main result of this paper, the radial scale factor and radial derivatives 
of other metric functions. In the previous section we discussed the solution t{a). Since we know 
by definition the exact one can solve numerically for a{t) using a simple Newton- Raphson 
algorithm, obtaining a{t) at machine accuracy level at hardly any computational cost. Therefore, 
in the following we assume {r, t\ as input parameters, and the function a(r, t) as a known function. 
We aim to express all solutions in terms of those quantities. 

The functions appearing in the metric are R{r,t) and R'ir.t). The time derivatives of these 
functions are relevant when one wants to solve geodesic equations in this metric, which is why we 
discuss them here as well. We have R{i\ i) = ra{r, t), such that 



R'{r,t) ^a{r,t) + ra'{r,t) 
R{r,t) ^ra{r,t) 
R'ir,t) =a{r,t) +ra{r,t) 

=d + ra'{r, t)H{r, t) + ra{r, t)H'{r, t) 



H'{r,t) 



1 



2H{r,t) 

1 / SttM^ r 



2H{r,t) \ 3 



-3 -3K(r) 
i(r,i) ^ 27ra3(r,t) 



a'{r,t) + 



(35) 
(36) 

(37) 
(38) 



Comparing Eq. ([571) to Eq. 



we see that we only need the solution for a'(r, t) in order to be 



able to calculate all relevant quantities. However, the term 2H\r t) ^^'^^ care to be taken at the 
transition from expansion to collapse, where H{r, t) = 0. We will show in the following that this 
limit is in fact regular, and that all metric functions remain properly defined throughout all the 
expansion and collapse history. 



5.1 Spatial derivative of the scale factor during expansion 

Since t is one of the orthogonal coordinates, we have drt = 0, even if we solve for t hy t 
t{a{r,t),r). Therefore, 

d d d 

— [t{a, r) - tBB{r)] ^a{r, t)—t{a, ^) + ^ [tia, r) - tBsir)] 



= - drtBsir) =a'(r,t)— ^ + |- [t(a, r) - tBsir)] 
a[r,t) or 

d 



a\r,t) = - d{r,t) 

The only non-trivial term in (j40[) is 

d d f 

M-[tia,r)-tBBir)]^- 



— [t{a, r) - tBB{r)] + drtBsir) 
or 



fa da 



Y' 



VXTYa^+Za? 
A r 

a 

X + Ya + Za^ 







da. 



(39) 
(40) 



(41) 

Y ~ 2n{r) and 

Z = -—n- This expression is again an elliptic integral. We spare the reader the detailed steps, but 



where we continue to use the notation of the previous section, A ~ a{r, t), X 

A 
3M2 



8lT 
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the general procedure is very much like in Eq (|28p . however one not only substitutes a 
and c ^ b = c — ^, but one also splits into partial fractions, to arrive at 



Y' 
2Zi 



Y' 



3 



§ + ^ a + a3 



da 



3Z2 (21222:3)- 
1 



1 



(22 - 21) (22 - 23) 
1 

(23 - 2l)(23 - 22) 



1 1 1 

{Zl ~ Z2){ZI ~ Z^i^^ \A Z2 A 

1 11 11 1 

A 23 ' ^ zi' A 22 
1 11 11 1 

Z3 



1 1 



1 



Rd 



A zi A Z2 A 



(42) 



One can take this equation one step further, using Eq. ([5]), which is R]j{x,y, z) + R]j{z,x,y) + 
R]j(y, z,x) = 3{xyz)~^. Hereby one eliminates one evaluation of the function R]j{x,y, z), and 
more importantly, it reveals the kind of singularity that is encountered when one of the arguments 
goes to zero. Choosing 21 to be the smallest positive real root, or as in the previous section 
U{r) = 21, we finally have the solution. 



Ma'{r,t) 
Y' 



A 



Z^iA- zi){A- Z2){A- Z3) 



Z [zi - 22)(2l - 23) 
I f 1 



A 



3Z2 (212223)2 
Y' i 



1 



(22 - 2l)(22 - 23) (21 - 22)(2l - 23) 
1 1 



3Zi (2122^3)5 V(^3 - Zl){z3 - 22) (21 - 22)(2l - 23) 
+MdrtBBir) 

with A = a(r, t), which can be recast as 

Ma'{r, t) = a(r, t)Q{r, t) + P(r, t), 

, , Y' A 
P{r,t) 



Rd 
Rd 



1 

1 



1 1 

23' A 
1 1 



1 1 1 

zi A 22 
11 1 

22 ' A 23 



Q{r,t) 
Y' 



Z [zi - Z2){ZI - 23) 

MdrtBB{r) 

1 



(43) 

(44) 
(45) 



+- 



3Z2 (212223): 
Y' I 



1 



(22 - 2l)(22 - 23) (21 - 22)(2l - 23) 
1 1 



2,Zi (2122^3)^ V(^3 - Zi)[zz - 22) (21 - 22)(2i - 23) 



Rd 
Rd 



1 


1 1 


1 


1 


1 


A " 


23' 








1 


1 1 


1 


1 


1 


1 ^ 


21 ' A 


Z2 




23 



(46) 



with as before A = a{r, t). 

If we look back at Eq. (|40| and Eq. (|42l) . we see that the overall factor d(r, t) in the expression 
for a'{r,t) multiplies inside dr [t{a,r) — tBB{r)], such that a'{r,t) is finite and non-zero for 
d{r,t) — )- 0, that is lim^(r ^^.^o a'(r, t) = P{r,t)/AI. 
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5.2 Spatial derivative of the scale factor during collapse 

During the collapsing phase, we have 



M[tiA)-tBB{r)] 



\fada 



\J{a- zi)(a - zi){a - Z3) 
^ ^Jada 
Z Jo ^(a - zi){a - Z2){a - 2:3) 



(47) 



We write explicitly the r-dependence in U (r) = (jlj , to point out that the transition from 
expansion to collapse is a priori an r-dcpcndcnt event jj Taking the derivative of this expression, 
we find 



Ma'{r,t)=\d\ 



dr 



21 W 



fa da 



XL 
T 



-/Z Jo ^/(a - zi){a - Z2){a - Z3) 
da + MdrtBBir) 



X + Ya + Za^ 



(48) 



where the absolute value |d| is to remind us that we take the positive root of \/{a — zi){a — Z2){a — Z3). 
When we realize that the r-dependence in zi{r) is entirely specified by Y{r) as that is the only 
r-dependent function in all integrals, such that drZi = Y' {r)dY zi, the first term becomes. 



dr 



fa da 



a/Z Jo yjia- zi){a - Z2){a - Z3) 

a.yiim^^-^^ ^ 

A^zi 



Z ^ [A - zi){A - Z2){A - Z'i) zi Jo \{a - zi){a - Z2){a - z:^) 



da 



(49) 



As zi satisfies X + Y zi + Z z\ ~ wc have dyzi 
Y + 3Zzf =Z lim Db (^ + ^B + B 



B^zi 



z 



Y+3Zzf 



fi-^- Next, we observe that 



=Z lim dB{B - Z2){B - Z3){B - zi) 

=Z lim [{B - Z3){B - zi) + (B - Z2)iB - zi) + (B - Z2)(S - Z3)] 

B— i.21 

= Z(Z1 - Z2)(Z1 - Z3). 



(50) 



so that dvzi = ^7 Y7 T- Using this relation, together with the solution to the same integral 

Z[Zi—Z2){Zi—Z^) O 7 o o 

in Section [5. 1[ and some more simple algebra, we arrive at 



<,ii(r, t)=d (2Q(r, tu) - Q{r, t)) + P{r, t), 



(51) 



where tu = t{a = Zi, r), the subscript 'coll' denotes that this expression is valid during a collapsing 
phase, Q{r, t) and P(r, t) are defined in Eqs. (|45l46p . and Q{r, tu) is evaluated by replacing A — s- zi 
in Eq. 

It should be understood now that with expressions ((44]) and (|5ip . a'(r, t) is finite and continuous 
at d 0. 



^Actually, we have U{r) = a{r,tij), where tu denotes the time at which d = 0. Since this time tjj is itself a 
function of r, one can show that U'{r) = a' {r, tu) + a(r, tij)drtij = drZi. 
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5.3 The spatial derivative of the Hubble parameter at a = 

During expansion, the expression for H'{r,t) in Eq (|38p is regular. At a — > 0, we can now insert 
Eq. [33]for a'{r,t), to find after some manipulations, 



3XY' 



a — zi 



* 2a3 aV 2x^a^3X + 2Yzi)\l {a- Z2){a- Z2)' 



(52) 



which is pefectly regular for the whole domain a(r, t) G {0, zi}. One extra minus sign appeared in 
the last expression, following from = —^/a — Z\ with a < z\ for all a. 

5.4 The spatial derivative of the Hubble parameter during collapse 

Combining Eq. (|5T|) and Eq. (|52|) . wc find during a collapsing phase, 

iXY' 



V2a^ aV 2V2^ai(3X + 2yzi) V (a-^2)(a-Z3) 

preserving continuity at d = 0. 



, (53) 



6 Example of application 

In Figure [I] we show an example of an application of the solutions presented in this work. The 
figure shows a comparison between two different density profiles, both parametrized by 



where K{, is defined in Eq. p7p . and with 

r 1 



1 

47r- 



1 

4ir^ 



1 



8 



r—aL 
L-aL 



-1 + 8^2 



r — aL 
L-aL 







cos ( 47r 



■ r~aL 
L-aL 



for r < aL 



for aL < r < 



for <r < L 

for r > L, 



(54) 



(55) 



which is the third order of the interpolating function W^„(x,a) which interpolates from 1 to in 
the interval a < a; < 1, while remaining C" everywhere. We introduce the interpolating function 
Wn{x,a) in Appendix ICl For < a < 1 this function is on r G [0, 00). Hence, all functions in 
the metric are everywhere, including at the center (r = 0) and at the matching to FLRW at 
r = L, guaranteeing a finite and continuous Ricci scalar everywhere. By construction all functions 
shown in the figure in fact have continuous first derivatives in r, even if by eye it may seem 
otherwise. 

The left column in Figure [T] shows quantities for a profile with a = 0, the right column shows 
the same quantities for a = |. In the top row we see the curvature profiles as a function of 
comoving radius, which are time independent. Both profiles share a same matter density at the 
centre where r = and at the outer radii r > L, however the difference in shape of profile leads to 
a different total amount of matter inside the over-density. For a = |, the over-density possesses 
a larger region with large closed curvature. The second row shows a time-dependent auxiliary 
radius. 



rFL-Rwir,t) 



S{r',t) 
a{oo,t)' 



dr\ 



(56) 
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Figure 1: Comparison of the two distinct overdensities. Curvature profile (top row) as a function 
of comoving radius, the auxihary FLRW radius rpLRW as a function of radius and time (second 
row), the matter density as a function of the auxihary radius rpLRW a-nd time (third row) and 
finally the time derivative of the metric function 5'(r, t) as a function of auxiliary radius rpLRW and 
time (bottom row). Both over-densities are matched to the same homogeneous ACDM- universe, 
at the same comoving radius r = 1 Mpc. The difference between the over-densities is in the 
value a in Eq. (|55p. determining the range in r over which the curvature profile falls back to the 
background value and thereby determining the total mass in the over-density. In all graphs, the 
time coordinate is represented by the colour of the curves, evolving from red to white to blue (from 
black to white to black in black-and-white print), indicating time varying from today (to = 13.3 
Gyr) to some moment in the past {to = 2.8 Gyr). Additionally, labels inside the graphs indicate 
times corresponding to different curves. 
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which keeps the length measure along this radius at a given time constant, ds^ = S(r,t)dr'^ = 
^''flrW' ^"-"^ constant in time. This radius can be interpreted loosely as what an FLRW 
observer would see. For these scenarios it turns out that roughly rpLRw — R'{r,t)/a{oo^t), which 
is its definition in Rcfs [Hs. 57|. 



In the third row we show the relative matter density, normalized to the matter density of the 
surrounding homogeneous cosmology. While at the centre the time evolution is the same for both 
profiles, the surrounding under-dense (but still closed curved) shell differs largely between the two 
cases. Since K{r) stays large and negative up to higher radius in the a = |-casc, compared to the 
a = 0-case, a larger range in r is present for which shells have a collapsing solution and actually 
experience the collapse. Therefore, the range in r in which the shells expand rapidly and become 
more and more under-dense, is smaller for the a = |-case. Hence, the resulting surrounding 
under-dcnsity is less dense and more emphasized in the a = |-case than in the a = 0-case. 

For a radial null geodesic, the trajectory is defined by, 

^^-S{r,t) (57) 
ar 

dz 

^ = (l + z)5(r,i). (58) 

This is why, for illustration, we plot S{r,t) in the fourth row in Figure [TJ 

Obviously, as the effect of the spherical collapse on the outer radii is more violent for the 
a = |-case, the red shift that a photon experiences by passing through that region is larger than 
in the a = 0-case. Even though this observation is not enough to draw conclusions about photon 
geodesies and collapsing structures, it illustrates how the solutions in this work can be used to 
further asses the importance of the initial distribution of matter in the line of sight, on distant 
observations. 

Note that all quantities remain perfectly smooth and continuous at the transition from expan- 
sion to collapse, which in the third and fourth row in Figure [1] occurs roughly where each curve 
crosses the level of the background, i.e. the level of the same curve at r > L. 

The solutions presented in this work allow for practically instantaneous calculation of the 
quantities presented in Figure [TJ We release a module, written in Fortran, which returns all 
metric functions and derivatives thereof as a function of time for given functions K(r) and tsBir), 
and given cosmological parameters. The module is released at 



http : / / web . physik . rwth-aachen . de/ download/ valkenburg/ ColLambda/ 



7 Conclusion 

We have presented an as complete as possible overview of the solutions to the Einstein equations 
governing the Lemaitre-Tolman-Bondi metric, including fully continuous solutions for collapsing 
over-densities surrounded by an expanding universe. The solutions are written in terms of Elliptic 
Integrals in Carlson's symmetric form, which allow for fast numerical evaluation of the solutions 
at machine accuracy level. The solutions to all metric functions involve the numerical inversion of 
one function, t{a), whose derivative is explicitly known a priori, therefore allowing for inversion 
at machine accuracy level while remaining sufficiently fast. 

We finished with a brief example of how these solutions can be applied. These solutions could 
improve the accuracy and speed of many analyses involving structure formation and inhomoge- 
neous cosmologies. 
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A Asymptotic expansions 



The solutions of the scenario where all three components, dust, curvature and the cosmological 
constant are non-zero and not asymptotically small, are given in the main body of this paper. 
For several purposes, a not unimportant one being numerical accuracy, asymptotic expansions are 
useful. By asymptotic expansions we mean the solutions for a given size of the scale factor, where 
one or more of the constituents contribute only marginally to the result. 

In the following we use a looser definition of fi^, where at each size of the scale factor fli denotes 
the fractional contribution to H{a). That is, for matter for example we re-define VLra{a,r) = 
and so on. 

We focus only on solutions with a Big Bang, and only solutions with a non-zero matter (dust) 
content. Then, given Eqs. ([TG]) and (|30)) . which we repeat here, 




=H (r, t) = — - 



(59) 



, , , , , 1 /" a/o da 

t{A)-tBB{r)=^l , ^ (60) 

M Jo §-. + 2K{r)a+ -^d^ 



we see that for any choice of K{r) and A. prior to some initial time the equation is dominated 
by matter. Hence the integral that gives t{a) always has a non- negligible contribution from the 
matter content, even for a so large that ^ g^J^j-y ■ Therefore we do not consider asymptotic 
expansions for r2,„ <C 1. 

Similarly, when expanding the integrand in Eq. (|60p for small omega's, it only makes sense to 
expand in ilk if ^ ^m, regardless of flA. To summarize, one can use expansions under the 
conditions, 

expand in fife when <e, (61) 

An 

Aa^ 

expand in fl\ when = <e, (62) 

8nAP + 6K{r)a 

guaranteeing that the expansion is valid throughout the whole integration from to a, for t{a) . 
When one uses a linear expansion, one should set e to e = ^/rj, with rj the desired accuracy, such 
that the error is C(e^) = 0{ri). In the case of numerical computation, one sets rj to the machine 
precision, which for double precision (64 bit floating point) is t] = 10~^^, such that e = 10~^. That 
is, in double precision one has approximately 16 significant digits. 

A.l fiA<l, fi„^0^1]fe 

For small A, the integral in Eq. ([BO]) becomes 



1^ Vdda A /"^ ai da / A^ \ 

M\t{A) -tBsir)] ^ , = ^/ -+0[-^] 

AkM r—r- — 2in . , _i / l3AK(r)\ ] 

^V3A.(0+4.-— smh ^(^^-Aijj 

_9 j VG^JAMt) {9A^n{rf - 2l7r^2^(r)2 + 70n^ AK{r) + 280n^) 




GAP'^'"^^ ' \ 108 ^ 3 AK{r) + An 

_g.3^.„.-.(;2f5)u„(^), 
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which reduces to the well-known result as A — )■ 0. The turning point from expansion to collapse 
lies at a = — g^^,.-, if this expression is positive. Care must be taken with the branch cut of the 



2 . For K 2 the minus signs cancel, and no care 



square root here, as for negative k, 
has to be taken. 

To obtain the derivative of time with respect to radius r in this case, one must first take the 
derivative of the full expression, and then expand, to arrive at. 



Mdr[t{A)-tBB{r)] = 



a'{r,t)M 
d{r,t) 



-TTsinh 




l3AK{r 



A K'{r 

M2V2K(r)5| 48(3AK(r) +47r)i 




(64) 



A.2 < 1, fi„ ^0 ^f^A 

For small K{r), the integral in Eq. becomes 



M[t{A)-tBB{r)] = 



I da 



Stt 
3 



3A/2" 



-Kir) 



8tt_ 
3 



A 53 



dd + O {nir f) 



(65) 



The first integral is in principle elliptic, and the full solution in Eq. ([31]) is applicable when one 
defines the correct roots Zi, however this special case is has a known solution in terms of sinh(a::) 
and can be found in the literature. The second integral is of course the integral that is solved in 
the main text for t'{a, r) for the special case where K{r) = 0. Hence, 



M[tiA)-tBB{r)] 



~3A" 



■ sinh 



2n{r) / A y i 



1 



(X2 - Xi){x2 - X3) 
1 

(X3 - Xi){xs, - X2) 



Rd 



Rj 



'a 

1 




1 



Rd 



[Xi - X2){xi - Xs) 
11 11 1 

xz^ A xi ' A X2 
11 11 1 
xi^ A X2' A 



1 1 

~2A 



1 1 



1 

Xl 



(66) 



where Xi are the three solutions of ^ 



j^x^ — 0. If one allows for a negative A, also this 
scenario can have a postive real Xi at which the transition from expansion to collapse occurs. Of 
course care has to be taken at the limit of a — > 0, identical to what is described in the main text 
concerning a'{r,t). 

Unfortimately, an expansion to obtain 9,. [t{A) — t_B_B(r)] in this asymptotic region is not trivial, 
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as, 



Mdr[t{A)-tBB{r)] = 

2K'{r) ( A 



a'{r,t)M 
a{r,t) 

i 



iAP J (xi 2:2 0:3) 2 



1 11 11 1 

{xi- X2){xi- xz)""^ \A X2' A xa'A xi 



-Rd 



1 11 11 1 



{X2 - Xi){x2 - X-i) \A X'i' A Xi' A X2 



1 11 11 1 

{xz- xi){x-i- X2)"'^^ \A xi' A X2' A X3 



-Rd 







f 


a 


Jo 


. 3 ^ 3M2 " . 



(67) 



where at this moment we do not know of a solution of the last integral in terms of symmetric 
elliptic integrals, so we leave that integral for future work. 



A.3 < 1, fife < 1, fi^ 7^ 

The simplest of the expansions is the scenario with small K(r) and A. The integral becomes. 



M[t{A)-tBB{r)] =\ — / d& 
87r J„ 



r 3 , ,.3 3 A ,7 

'a Klrja'^ = — 

Stt ^ ^ 87r6M2 



with the spatial derivative 

Mdr[t{A) - tBB{r) = 
a'{r,t)M 



a{r,t) 



A 



Stt/ 11M2 



(69) 



B Solution in presence of radiation 

Writing A = a(r, t), the general solution to t for the Fricdmann equations in presence of radiation, 
matter, curvature and a cosmological constant is given by the integral. 



I da 



H,[t{A)-tBB{r)] 



(70) 



which is the same as Eq. (p5)) for n = 4. 

We split the integral in two parts, where it is most convenient to take zi to be the smallest 
positive real root, if any root is positive and real (otherwise any root will do). 



[t{A) tBB] / 



{a — Zi + Zi) da 

° \/nl=i(«-^™) 

{a — zi) da 



A Jo 



(71) 
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In the second integral we can substitute a 



1 _ i_ 

a A ' 



da 



1 



111=1 (a - Zra) ^ \l nl=2 
1 



db 



Rf{Vi2, V23). 



(72) 



^here V,, = yjTTyTTX^ y^i - ^ with {^,J, k,l} ^ {I, 2, 3, 4} {e.g. if j = 2, 3, 

then fc,Z — 1,4). The last equality in Eq. ([7^ is proven in Ref. [t^, and the equality is invariant 
under the ordering of the roots z„,, that is, invariant under the choice for Z4. 

The first integral in Eq. ([71]) takes another road. Following Ref. [t^I we make the change of 

Ab 

+1 



variables a — )• 



1 



A Jo 



(a — Zl) da 



1 



A Jo 



Ab 
b+1 



Zl 



{{A - zi)b - Zl) db 



^aJo ib+l)Utn=lVi^-^m)b~Z, 



A - Zl 



^\Iltn=2i^-^rn)Jo (6+1)01=1 



A—zi 



db 



A-z 



1 



A~ Zl 



In the last line we use the definitions. 



A-z^ 



A- 



Zk 



Zl 



A - zkV A~ Zl 



with {z,j,fc,Z} = {1,2,3,4}, 



2 TT/2 ^ik^ili^j 1) 



(73) 

(74) 
(75) 
(76) 

(77) 

(78) 

In Eq. ((73)) we wrote the intermediate step in the second line explicitly, since from there it is 
straightforward to see that for real positive zi , we have 

'■"^ (a — Zl) da 



a =1 



'J Zl - 1 

Z2y/zEy/Z4 



lim 



r =-W^io- 

^1 



1 



2^ 



1 



Z2 



Z3 



Z4 



Zl - Zl Zl - Z-i Zl - Z4 

(79) 



The second integral in Eq. (j7ip is continuous in this limit. 
Alltogether this gives the final result. 



i/* \t{A) - tBB\ = 



1 



A - Zl 



"AVm=2(^-2m) 

1 xAT 



''i{zi-\) ^' ^^'3' + 



llm=2 -m 



(21 - 1) 
i?F(V\|,Fiiy23)- 



(80) 
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B.l Roots Zi 

In order to write down the roots {zj}, let us use the fohowmg definitions, 



X 



K 



L 



Y 



{-72XZ + 27Y^ -f 
2Z 2i{12X + Z'^ 



3K 



2Z^f - 4 {12X + Z'^f - 72XZ + 271"^ 



K 

233' 



2Z' 



(81) 
(82) 
(83) 



The roots are 



Z2 



Z4 




6Z 


- L- 


2Y 

7i' 


(84) 


6Z 


- L- 


2Y 

7l' 


(85) 


6Z 


-L + 


2Y 


(86) 


6Z 


-L + 


2Y 

7l- 


(87) 



C The interpolating function Wn{x, a) 

We define the function Wn{x,a) as, 
Wnix,a)= 1 



for n > 2. 



W2{x,a)= i + isin TT i-f^ 



W„+i{x,a) = 
Wn{x,a) = 



/o''° W„{\l-2x'\.0)dx' 
W,A\l~2x"\,0)dx" 



for X < a 



for a < X < 1 



for 1 < 



which is C" everywhere. Note the absolute value inside the integral, which takes care of mirroring 
two images of the function that interpolates from one to zero, in order to have a function that goes 
from zero to one to zero. Also, inside the integral the functions are evaluated for a = 0, because 
the integration limits are always < {x',x"} < 1. 

One can construct a similar function taking polynomials in stead of sinusoids, by starting with 
Wo{x,a) = (1 — x)/{l — a) for a < a; < 1. 

We constructed this function by starting with the function f{x) = 1 + sin(a;), which equals 
zero and has a zero derivate at x = — 7r/2 and at x = 3tt/2. Integrating this function over x then 
leads to a function which has zero first and second derivatives at x = —tt/2 and at x = 37r/2. 
Next one matches this function too a mirror image, normalizes it to zero at both end points of 
the domain and integrates again, such that one ends up with a function that has zero first, second 
and third derivatives at the domain borders. One can continue this scheme forever, as explicated 
in Eq. 



D Usage of the numerical module ColLambda 

The module ColLambda is written in Fortran90, and can be downloaded from 

http : // web . phys ik . rwth- aachen . de / download/ valkenburg/ ColLambda/^ Compilation instruc- 
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tions can be found at that website. The module can be included by the statement, 
use LLTB 

and any metric function and a number of derivatives is then obtained by the function call, 

call lltb_f unctions (HO_inf , Lambda, kofr, dkdr, tbbofr, dtbbdr, Mtilde, & 
r, t, & 

Rltb, Rp, Rpd, Rpdd, S, Sd, Sdd, a, ap, apd, add, apdd, H, Hp, tturn) 

where all items in the first and second line are mandatory and are intent (in) , and all items in the 
third line are optional and intent (out). The normalization parameters and their corresponding 
parameters in this paper are 

HO_inf ^H^=H^{n,to) 
Lambda — A 
Mtilde = 

and the functions 

kofr 

dkdr 
tbbofr 
dtbbdr 

must each be a pure function of r, declared as: 

function kofr(r) 

real (8) : : kofr 

real (8), intent (in) :: r 
end function kofr 

As a consequence, any other parameters on which k{r) may depend, such as the maximum size L, 
maximum curvature Kmax, or anything else, must be global variables which the function kofr(r) 
can access without receiving them as arguments. The values that the subroutine returns have the 
following notation: 



Rltb 


= i?(r,t) 


Rp 


= R'{r,t) 


Rpd 


= R'{r,t) 


Rpdd 


= R'{r,t) 


S 


= 5(r,t) 


Sd 


= S{r,t) 


Sdd 


= S{r,t) 


a 


= a{r,t) 


ap 


= a'{r,t) 


apd 


= a'{r,t) 


add 


= a{r,t) 


apdd 


= -d'{r,t) 


H 


^H{r,t) 


Hp 


= H'{r,t) 


tturn 


= tu{r) 







where tu{r) denotes the time at which a{r^tu) ~ 0, if it exists. If it exists, the local singularity is 
reached at t = tBsir) + itjj, that is, a{r,tBB{r)) ~ a{r,tBB{r) + 2tij) = 0. If this time does not 
exist, i.e. when the solution is ever expanding, tturn will be set to lO'^". 

An example call, where the user has set the normalization variables and has defined the neces- 
sary functions K(r), i_BB(^) and their first derivatives, to set myRltb to R{r,t), myS to S{r,t) and 
myt to tij{r), would look like: 

call lltb_f unctions (HO_inf , Lambda, kofr, dkdr, tbbofr, dtbbdr, Mtilde, & 
r, t, & 

Rltb=myRltb, S=myS ,tturn=myt) 



= K(r) 
= K'{r) 
^tBB{r) 
= 4bW 
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